home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / bilinear.pro < prev    next >
Text File  |  1997-07-08  |  4KB  |  127 lines
  1. ; $Id: bilinear.pro,v 1.3 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1993-1997, Research Systems, Inc.  All rights reserved.
  4. ;    Unauthorized reproduction prohibited.
  5.  
  6. FUNCTION BILINEAR,P,IX,JY
  7. ;+
  8. ; NAME:
  9. ;    BILINEAR
  10. ;
  11. ; PURPOSE:
  12. ;    Bilinearly interpolate a set of reference points.
  13. ;
  14. ; CALLING SEQUENCE:
  15. ;    Result = BILINEAR(P, IX, JY)
  16. ;
  17. ; INPUTS:                 
  18. ;    P:  A two-dimensional data array.
  19. ;
  20. ;    IX and JY:  The "virtual subscripts" of P to look up values
  21. ;      for the output.
  22. ;  
  23. ;    IX can be one of two types:
  24. ;         1)    A one-dimensional, floating-point array of subscripts to look
  25. ;        up in P.  The same set of subscripts is used for all rows in
  26. ;        the output array.
  27. ;         2)    A two-dimensional, floating-point array that contains both 
  28. ;        "x-axis" and "y-axis" subscripts specified for all points in
  29. ;        the output array.
  30. ;
  31. ;    In either case, IX must satisfy the expression,
  32. ;            0 <= MIN(IX) < N0  and 0 < MAX(IX) <= N0
  33. ;    where N0 is the total number of subscripts in the first dimension
  34. ;    of P.
  35. ;
  36. ;    JY can be one of two types:
  37. ;         1) A one-dimensional, floating-point array of subscripts to look
  38. ;        up in P.  The same set of subscripts is used for all rows in
  39. ;        the output array.
  40. ;         2) A two-dimensional, floating-point array that contains both
  41. ;               "x-axis" and "y-axis" subscripts specified for all points in
  42. ;               the output array.
  43. ;
  44. ;        In either case JY must satisfy the expression,
  45. ;            0 <= MIN(JY) < M0  and 0 < MAX(JY) <= M0
  46. ;        where M0 is the total number of subscripts in the second dimension
  47. ;        of P.
  48. ;
  49. ;      It is better to use two-dimensional arrays for IX and JY when calling
  50. ;      BILINEAR because the algorithm is somewhat faster.  If IX and JY are 
  51. ;      one-dimensional, they are converted to two-dimensional arrays on
  52. ;      return from the function.  The new IX and JY can be re-used on 
  53. ;    subsequent calls to take advantage of the faster, 2D algorithm.  The 
  54. ;    2D array P is unchanged upon return.
  55. ;
  56. ; OUTPUT:
  57. ;    The two-dimensional, floating-point, interpolated array.  
  58. ;
  59. ; SIDE EFFECTS:
  60. ;    This function can take a long time to execute.
  61. ;
  62. ; RESTRICTIONS:
  63. ;    None.
  64. ;
  65. ; EXAMPLE:
  66. ;    Suppose P = FLTARR(3,3), IX = [.1, .2], and JY = [.6, 2.1] then
  67. ;    the result of the command:
  68. ;        Z = BILINEAR(P, IX, JY)
  69. ;    Z(0,0) will be returned as though it where equal to P(.1,.6) 
  70. ;    interpolated from the nearest neighbors at P(0,0), P(1,0), P(1,1)
  71. ;    and P(0,1).
  72. ;
  73. ; PROCEDURE:
  74. ;    Uses bilinear interpolation algorithm to evaluate each element
  75. ;    in the result  at virtual coordinates contained in IX and JY with 
  76. ;    the data in P.                                                          
  77. ;
  78. ; REVISION HISTORY:
  79. ;       Nov. 1985  Written by L. Kramer (U. of Maryland/U. Res. Found.)
  80. ;    Aug. 1990  TJA simple bug fix, contributed by Marion Legg of NASA Ames
  81. ;    Sep. 1992  DMS, Scrapped the interpolat part and now use INTERPOLATE
  82. ;-
  83.     ON_ERROR,2              ;Return to caller if an error occurs    
  84.     IF((N_ELEMENTS(IX) EQ 0) AND (N_ELEMENTS(JY) EQ 0)) THEN BEGIN
  85.       I=FIX(IX) & J=FIX(JY) & IP=I+1 & JP=J+1
  86.       DX=IX-FLOAT(I) & DY=JY-FLOAT(J)
  87.       DX1=(1.-DX) & DY1=(1.-DY) 
  88.       RETURN,( P[I,J]*DX1*DY1 + P[I,JP]*DX1*DY $
  89.           + P[IP,J]*DX*DY1 + P[IP,JP]*DX*DY)
  90.     ENDIF
  91.  
  92.     A=SIZE(IX)  & B=SIZE(JY)
  93.     NX=A[1]
  94.     IF(B[0] EQ 1) THEN BEGIN
  95.       NY=B[1]
  96.     ENDIF ELSE BEGIN
  97.       NY=B[2]
  98.     ENDELSE                                                   
  99.         IF(A[0] EQ 1) THEN BEGIN
  100.       TEMP=IX
  101.       IX=FLTARR(NX,NY)
  102.       FOR I=0,NY-1 DO IX[0,I]=TEMP
  103.     ENDIF
  104.     IF(B[0] EQ 1) THEN BEGIN
  105.       TEMP=JY
  106.       JY=FLTARR(NY,NX)
  107.       FOR I=0,NX-1 DO JY[0,I]=TEMP
  108.       JY=TRANSPOSE(JY)
  109.     ENDIF
  110.     return, interpolate(p, ix, jy)    ;Use new interpolate function
  111. ;    I=FIX(IX) & J=FIX(JY)
  112. ;    IP=I+1   &  JP=J+1
  113. ;    DX=IX-I & DY=JY-J
  114. ;    DX1=1.-DX & DY1=1.-DY
  115. ;    Z=FLTARR(N_ELEMENTS(I(*,0)),N_ELEMENTS(J(0,*)))
  116. ;    NUMX=N_ELEMENTS(I)
  117. ;    PZ=FLTARR(N_ELEMENTS(P(*,0))+1,N_ELEMENTS(P(0,*))+1)
  118. ;    PZ(0,0)=P(0:*,0:*)
  119. ;    FOR N=0L,NUMX-1 DO BEGIN
  120. ;      Z(N)=  PZ(I(N), J(N)) *DX1(N)*DY1(N)    $
  121. ;    +        PZ(I(N),JP(N)) *DX1(N)*DY(N)    $
  122. ;    +        PZ(IP(N),J(N)) *DX(N) *DY1(N)    $
  123. ;    +        PZ(IP(N),JP(N))*DX(N) *DY(N)
  124. ;    ENDFOR 
  125. ;    RETURN,Z
  126. END            
  127.